clear all; clc; close all;
%Updated: 21/09/2018


cd 'baseline'

%v00b    =load('v.txt');
ctt00b  =load('gct.txt');
h00b    =load('gh.txt');
gn00b   =load('ggn.txt');
def00b  =load('gdef1.txt');
ibp00b   =load('gibp.txt');
vcon00b = load('gvcon.txt'); %v00b(:,2);
vaut00b = load('gvaut.txt'); %v00b(:,3);

q_paid00b =load('gq.txt');
qpaid00b  = q_paid00b(:,1);
qnodef00b = q_paid00b(:,1);

cd ..
cd 'commit/period1'

v00    =load('v.txt');
ctt00  =load('ctt.txt');
h00    =load('x.txt');
gn00   =load('gn.txt');
def00  =load('default.txt');
bp00   =load('b_next.txt');
vcon00 = v00(:,2);
vaut00 = v00(:,3);

cttaut00  =load('cttaut.txt');
haut00    =load('haut.txt');
gnaut00   =load('gnaut.txt');

q_paid00 =load('q_paid.txt');
qpaid00  = q_paid00(:,1);
%qnodef00 = q_paid00(:,2);


cd ..
cd 'period2'

v11    =load('v.txt');
ctt11  =load('ctt.txt');
h11    =load('x.txt');
gn11   =load('gn.txt');
def11  =load('default.txt');
bp11   =load('b_next.txt');
vcon11 = v11(:,2);
vaut11 = v11(:,3);

cttaut11  =load('cttaut.txt');
haut11    =load('haut.txt');
gnaut11   =load('gnaut.txt');

q_paid11 =load('q_paid.txt');
qpaid11  = q_paid11(:,1);
%qnodef11 = q_paid11(:,2);

%cd ..
%cd ..
b       = load('bgrid2.txt');
a       = load('agrid2.txt');

nb=3501; %size(b,1);
na=501; %size(vcon00,1)/nb;

omegac = 0.30;
mu     = 1;
alpha  = 0.75;
r      = 0.02;
lambda =0.184;

for i=1:na;
    haut1(i)     =haut11(i);
    gnaut1(i)    =gnaut11(i);
    cnaut1(i)    =haut1(i).^alpha-gnaut1(i);
            
    haut0(i)     =haut00(i);
    gnaut0(i)    =gnaut00(i);
    cnaut0(i)    =haut0(i).^alpha-gnaut0(i);
 
	vaut0b(i,:)  =vaut00b(i);
                
    for j=1:nb;
        def1(i,j)   =def11((j-1)*na +i);
        if (def1(i,j)<111);
            vcon1(i,j)  =vcon11((j-1)*na +i);
            vaut1(i,j)  =vaut11((j-1)*na +i);
            qpaid1(i,j) =qpaid11((j-1)*na +i);
            qnodef1(i,j)=0; %qnodef11((j-1)*na +i);
            r1(i,j)     =100*((1+qpaid1(i,j).^(-1)-lambda)./(1+r)-1);

            h1(i,j)     =h11((j-1)*na +i);
            gn1(i,j)    =gn11((j-1)*na +i);
            cn1(i,j)    =h1(i,j).^alpha-gn1(i,j);
            ct1(i,j)    =ctt11((j-1)*na +i);
            pn1(i,j)    =(1-omegac)/omegac*(ct1(i,j)/cn1(i,j)).^(1+mu);
            bp1(i,j)    =b(bp11((j-1)*na +i));
        else;
            vcon1(i,j)  =NaN;
            vaut1(i,j)  =vaut11((j-1)*na +i);
            qpaid1(i,j) =NaN;
            qnodef1(i,j)=NaN;
            r1(i,j)     =NaN;
            h1(i,j)     =NaN;
            gn1(i,j)    =NaN;
            cn1(i,j)    =NaN;
            ct1(i,j)    =NaN;
            pn1(i,j)    =NaN;
            bp1(i,j)    =NaN;
        end;
                
        def0(i,j)   =def00((j-1)*na +i);
        if (def0(i,j)<111);            
            vcon0(i,j)  =vcon00((j-1)*na +i);
            vaut0(i,j)  =vaut00((j-1)*na +i);
            qpaid0(i,j) =qpaid00((j-1)*na +i);
            qnodef0(i,j)=0; %qnodef00((j-1)*na +i);
            r0(i,j)     =100*((1+qpaid0(i,j).^(-1)-lambda)./(1+r)-1);
            h0(i,j)     =h00((j-1)*na +i);
            gn0(i,j)    =gn00((j-1)*na +i);
            cn0(i,j)    =h0(i,j).^alpha-gn0(i,j);
            ct0(i,j)    =ctt00((j-1)*na +i);
            pn0(i,j)    =(1-omegac)/omegac*(ct0(i,j)/cn0(i,j)).^(1+mu);
            bp0(i,j)    =b(bp00((j-1)*na +i));
        else
            vcon0(i,j)  =NaN;
            vaut0(i,j)  =vaut00((j-1)*na +i);
            qpaid0(i,j) =NaN;
            qnodef0(i,j)=NaN;
            r0(i,j)     =NaN;
            h0(i,j)     =NaN;
            gn0(i,j)    =NaN;
            cn0(i,j)    =NaN;
            ct0(i,j)    =NaN;
            pn0(i,j)    =NaN;
            bp0(i,j)    =NaN;
        end;
           end;
end;


%%

for i=1:na,
    for j=1:nb;
        
        def0b(i,j)   =def00b((i-1)*nb +j);
        if (def0(i,j)==1);            
            vcon0b(i,j)  =vcon00b((i-1)*nb +j);
            qpaid0b(i,j) =qpaid00b((i-1)*nb +j);
            qnodef0b(i,j)=qnodef00b((i-1)*nb +j);
            r0b(i,j)     =100*((1+qnodef0b(i,j).^(-1)-lambda)./(1+r)-1);
            h0b(i,j)     =h00b((i-1)*nb +j);
            gn0b(i,j)    =gn00b((i-1)*nb +j);
            cn0b(i,j)    =h0b(i,j).^alpha-gn0b(i,j);
            ct0b(i,j)    =ctt00b((i-1)*nb +j);
            pn0b(i,j)    =(1-omegac)/omegac*(ct0b(i,j)/cn0b(i,j)).^(1+mu);
            bp0b(i,j)    =b(ibp00b((i-1)*nb +j));
        else
            vcon0b(i,j)  =NaN;
%            vaut0b(i,j)  =vaut00b((i-1)*nb +j);
            qpaid0b(i,j) =NaN;
            qnodef0b(i,j)=NaN;
            r0b(i,j)     =NaN;
            h0b(i,j)     =NaN;
            gn0b(i,j)    =NaN;
            cn0b(i,j)    =NaN;
            ct0b(i,j)    =NaN;
            pn0b(i,j)    =NaN;
            bp0b(i,j)    =NaN;
        end;
 
    end;
end;

%b=-b/(lambda+r);
%bp0=-bp0./(lambda+r);
%bp1=-bp1./(lambda+r);


lowa = round(na/4)+3;   %uncond mean-1std dev
meda = round(na/2);
hgha = round(3*na/4)-3; %uncond mean -1std dev

figure;
subplot(1,2,1);contourf(b(1:nb),a(1:na),def0(:,:),20);
xlabel('b');
ylabel('y^T');
subplot(1,2,2);contourf(b(1:nb),a(1:na),def1(:,:),20);
xlabel('b');
ylabel('y^T');

% ig=10;
% 
% figure;
% subplot(2,2,1);
% plot(b(1:nb),vcon_22(1:5:na,:,ig));
% xlabel('b'); title('g^N period 22') ;
% subplot(2,2,2);
% plot(b(1:nb),h_22(1:5:na,:,ig));
% xlabel('b');title('h period 22') ;
% subplot(2,2,3);
% plot(b(1:nb),ct_22(1:5:na,:,ig));
% xlabel('b');title('c^T period 22') ;
% subplot(2,2,4);
% plot(b(1:nb),bp_22(1:5:na,:,ig));
% xlabel('b');title('bprime period 22') ;
% dbstop;


% figure;
% subplot(3,2,1);
% plot(b(1:nb),vcon0(1:5:na,:)-vcon0b(1:5:na,:),'LineWidth',2);axis tight;
% xlabel('debt'); title('v^R period 1') ;
% subplot(3,2,2);
% plot(b(1:nb),h0(1:5:na,:)-h0b(1:5:na,:),'LineWidth',2);axis tight;
% xlabel('debt');title('h period 1') ;
% subplot(3,2,3);
% plot(b(1:nb),r0(1:5:na,:)-r0b(1:5:na,:),'LineWidth',2);axis tight;
% xlabel('debt');title('spreads period 1') ;
% subplot(3,2,4);
% plot(b(1:nb),bp0(1:5:na,:)-bp0b(1:5:na,:),'LineWidth',2);axis tight;
% xlabel('debt');title('bprime period 1') ;
% subplot(3,2,5);
% plot(b(1:nb),ct0(1:5:na,:)-ct0b(1:5:na,:),'LineWidth',2);axis tight;
% xlabel('debt');title('c^T period 1') ;
% subplot(3,2,6);
% plot(b(1:nb),gn0(1:5:na,:)-gn0b(1:5:na,:),'LineWidth',2);axis tight;
% xlabel('debt');title('g^N period 1') ;
%%
xx=1;
da=101;
db=1:50:nb;


figure;
subplot(3,2,1);
plot(b(db),vcon0(da,db)-xx*vcon0b(da,db),'LineWidth',2);axis tight;
xlabel('debt'); title('v^R period 1') ;
subplot(3,2,2);
plot(b(db),h0(da,db)-xx*h0b(da,db),'LineWidth',2);axis tight;
xlabel('debt');title('h period 1') ;
subplot(3,2,3);
plot(b(db),r0(da,db)-xx*r0b(da,db),'LineWidth',2);axis tight;
xlabel('debt');title('spreads period 1') ;
subplot(3,2,4);
plot(b(db),bp0(da,db)-xx*bp0b(da,db),'LineWidth',2);axis tight;
xlabel('debt');title('bprime period 1') ;
subplot(3,2,5);
plot(b(db),ct0(da,db)-xx*ct0b(da,db),'LineWidth',2);axis tight;
xlabel('debt');title('c^T period 1') ;
subplot(3,2,6);
plot(b(db),gn0(da,db)-xx*gn0b(da,db),'LineWidth',2);axis tight;
xlabel('debt');title('g^N period 1') ;

% figure;
% subplot(3,2,1);
% plot(b(1:nb),vcon1(1:5:na,:)-vcon0b(1:5:na,:),'LineWidth',2);axis tight;
% xlabel('debt'); title('v^R period 2') ;
% subplot(3,2,2);
% plot(b(1:nb),h1(1:5:na,:)-h0b(1:5:na,:),'LineWidth',2);axis tight;
% xlabel('debt');title('h period 2') ;
% subplot(3,2,3);
% plot(b(1:nb),r1(1:5:na,:)-r0b(1:5:na,:),'LineWidth',2);axis tight;
% xlabel('debt');title('spreads period 2') ;
% subplot(3,2,4);
% plot(b(1:nb),bp1(1:5:na,:)-bp0b(1:5:na,:),'LineWidth',2);axis tight;
% xlabel('debt');title('bprime period 2') ;
% subplot(3,2,5);
% plot(b(1:nb),ct1(1:5:na,:)-ct0b(1:5:na,:),'LineWidth',2);axis tight;
% xlabel('debt');title('c^T period 2') ;
% subplot(3,2,6);
% plot(b(1:nb),gn1(1:5:na,:)-gn0b(1:5:na,:),'LineWidth',2);axis tight;
% xlabel('debt');title('g^N period 2') ;


figure;
subplot(3,2,1);
plot(b(db),vcon1(da,db)-xx*vcon0b(da,db),'LineWidth',2);axis tight;
xlabel('debt'); title('v^R period 2') ;
subplot(3,2,2);
plot(b(db),h1(da,db)-xx*h0b(da,db),'LineWidth',2);axis tight;
xlabel('debt');title('h period 2') ;
subplot(3,2,3);
plot(b(db),r1(da,db)-xx*r0b(da,db),'LineWidth',2);axis tight;
xlabel('debt');title('spreads period 2') ;
subplot(3,2,4);
plot(b(db),bp1(da,db)-xx*bp0b(da,db),'LineWidth',2);axis tight;
xlabel('debt');title('bprime period 2') ;
subplot(3,2,5);
plot(b(db),ct1(da,db)-xx*ct0b(da,db),'LineWidth',2);axis tight;
xlabel('debt');title('c^T period 2') ;
subplot(3,2,6);
plot(b(db),gn1(da,db)-xx*gn0b(da,db),'LineWidth',2);axis tight;
xlabel('debt');title('g^N period 2') ;
dbstop;


figure;
subplot(2,2,1);
plot(b(1:db:nb),vcon(da/2,1:db:nb));
xlabel('debt');title('vcon low cT') ;
subplot(2,2,2);
plot(b(1:db:nb),vcon(na/2:da:na,1:db:nb));
xlabel('debt');title('vcon high cT');
%plot(b(1:nb),vcon(lowa,:),'-r',b(1:nb),vcon(meda,:),':r',b(1:nb),...
%     vaut(lowa,:),'-b',b(1:nb),vaut(meda,:),':b','LineWidth',1.5);
%xlabel('debt');title('value functions ');
% axis([b(1) b(nb) -35 -15]);
subplot(2,2,3);
plot(b(1:nb),h(1:2:na/2,:));
xlabel('debt');title('h low cT') ;
subplot(2,2,4);
plot(b(1:nb),h(na/2:2:na,:));
xlabel('b');title('h high cT') ;



figure;
subplot(2,2,1);
plot(b(1:nb),gn(1:2:na/2,:));
xlabel('b');title('gn low cT') ;
subplot(2,2,2);
plot(b(1:nb),gn(na/2:2:na,:));
xlabel('b');title('gn high cT');
%plot(b(1:nb),vcon(lowa,:),'-r',b(1:nb),vcon(meda,:),':r',b(1:nb),...
%     vaut(lowa,:),'-b',b(1:nb),vaut(meda,:),':b','LineWidth',1.5);
%xlabel('b');title('value functions ');
% axis([b(1) b(nb) -35 -15]);
subplot(2,2,3);
plot(b(1:nb),ct(1:2:na/2,:));
xlabel('b');title('ct low cT') ;
subplot(2,2,4);
plot(b(1:nb),ct(na/2:2:na,:));
xlabel('b');title('ct high cT') ;
dbstop;
dbstop;

figure;
subplot(2,1,1);
plot(a(1:na),vcon(:,1:2:nb/2),a(1:na),vaut(:,1),':');
xlabel('a');
subplot(2,1,2);
plot(a(1:na),vcon(:,nb/2:2:nb),a(1:na),vaut(:,1),':');
xlabel('a');
%plot(b(1:nb),vcon(lowa,:),'-r',b(1:nb),vcon(meda,:),':r',b(1:nb),...
%     vaut(lowa,:),'-b',b(1:nb),vaut(meda,:),':b','LineWidth',1.5);
%xlabel('b');title('value functions ');
% axis([b(1) b(nb) -35 -15]);


figure;
contourf(b(1:nb),a(1:na),def(:,:),20);
xlabel('b');
ylabel('a');



figure;
subplot(2,2,1);
plot(b(1:nb),bp(1:na/2,:));
xlabel('b');
subplot(2,2,2);
plot(b(1:nb),bp(na/2:na,:));
xlabel('b');
subplot(2,2,3);
plot(b(1:nb),q(1:na/2,:));
xlabel('b');
subplot(2,2,4);
plot(b(1:nb),q(na/2:na,:));
xlabel('b');
%plot(b(1:nb),bp(lowa,:),'-r',b(1:nb),bp(meda,:),':r',b(1:nb),bp(hgha,:),'--r','LineWidth',1.5);
xlabel('b');
title('debt policy');



% a   =load('a.txt');
% b   =load('b.txt');
% 
% 
% p0  =load('gp.txt');
% h0  =load('gh.txt');
% ibp0=load('gibp.txt');
% ct0 =load('gct.txt');
% m0  =load('gm.txt');
% gn0 =load('ggn.txt');
% tau0=load('gtau.txt');
% w0  =load('gw.txt');
% def0=load('gdef1.txt');
% q0  =load('gq.txt');
% 
% paut0  =load('gpaut.txt');
% haut0  =load('ghaut.txt');
% ctaut0 =load('gctaut.txt');
% maut0  =load('gmaut.txt');
% gnaut0 =load('ggnaut.txt');
% tauaut0=load('gtauaut.txt');
% waut0  =load('gwaut.txt');
% 
% na=rows(a);
% nk=1;
% ns=na/nk;
% nb=rows(b);
% 
% alpha = 0.63;
% kappa  =0.133;
% pm     = 1;
% gamma  = 0.43;
% omegac = 0.3;
% yloss  = 0.98;
% 
% atrans = ones(nb,1)*a';
% mstar = (pm./(a.*gamma)).^(1/(gamma-1));
% 
% for i=1:na;
%     for j=1:nb;
%         vcon(i,j)  =vcon0((i-1)*nb + j);
%         
%         p(i,j)  =p0((i-1)*nb + j);
%         h(i,j)  =h0((i-1)*nb + j);
%         ibp(i,j)=ibp0((i-1)*nb + j);
%         ct(i,j) =ct0((i-1)*nb + j);
%         m(i,j)  =m0((i-1)*nb + j);
%         gn(i,j) =gn0((i-1)*nb + j);       
%         tau(i,j)=tau0((i-1)*nb + j);       
%         w(i,j)  =w0((i-1)*nb + j);   
%         def(i,j)=def0((i-1)*nb + j);   
%         q(i,j)  =q0((i-1)*nb + j);   
%         rec(i,j)=w(i,j)*h(i,j)*tau(i,j);
%     end;
% end;
% 
% cn =h.^alpha-gn;
% const = (ct.^omegac).*(cn.^(1d0-omegac));
% bp = b(ibp);
% 
% for i=1:na;
%     vaut(i)  =vaut0(i);
%     
%     paut(i)  =paut0(i);
%     haut(i)  =haut0(i);
%     ctaut(i) =ctaut0(i);
%     maut(i)  =maut0(i);
%     gnaut(i) =gnaut0(i);       
%     tauaut(i)=tauaut0(i);       
%     waut(i)  =waut0(i);   
% end;
% 
% 
% cnaut =haut.^alpha-gnaut;
% constaut = (ctaut.^omegac).*(cnaut.^(1d0-omegac));
% 
% for i=1:na;
%     temp1 = find(ct(i,:) > 0);
%     temp2 = find(cn(i,:) > 0);
%     temp3 = find(p(i,:) > 0);
% 
%     iminct(i) = temp1(1);
%     imincn(i) = temp2(1);
%     iminp(i) = temp3(1);
%     
%     imin(i)  = max(iminct(i),max(imincn(i),iminp(i)));
%     
%     v(i,1:imin(i)-1)  = NaN; 
%     p(i,1:imin(i)-1)  = NaN; 
%     bp(i,1:imin(i)-1) = NaN; 
%     ibp(i,1:imin(i)-1)= NaN; 
%     ct(i,1:imin(i)-1) = NaN; 
%     cn(i,1:imin(i)-1) = NaN; 
%     const(i,1:imin(i)-1) = NaN;     
%     gn(i,1:imin(i)-1) = NaN; 
%     h(i,1:imin(i)-1)  = NaN; 
%     m(i,1:imin(i)-1)  = NaN; 
%     tau(i,1:imin(i)-1)= NaN; 
% 	rec(i,1:imin(i)-1)= NaN;
% %    def(i,1:imin(i)-1)= NaN; 
% 	q(i,1:imin(i)-1)  = NaN;
% 
% %     temp1 = find(ctaut(i,:) > 0);
% %     temp2 = find(cnaut(i,:) > 0);
% %     temp3 = find(paut(i,:) > 0);
% % 
% %     iminctaut(i) = temp1(1);
% %     imincnaut(i) = temp2(1);
% %     iminpaut(i)  = temp3(1);
% % 
% %     imin(i)  = max(iminctaut(i),max(imincnaut(i),iminpaut(i)));
% % 
% %     vaut(i,1:imin(i)-1)  = NaN; 
% %     paut(i,1:imin(i)-1)  = NaN; 
% %     ctaut(i,1:imin(i)-1) = NaN; 
% %     cnaut(i,1:imin(i)-1) = NaN; 
% %     constaut(i,1:imin(i)-1) = NaN;     
% %     gnaut(i,1:imin(i)-1) = NaN; 
% %     haut(i,1:imin(i)-1)  = NaN; 
% %     maut(i,1:imin(i)-1)  = NaN; 
% %     tauaut(i,1:imin(i)-1)= NaN; 
% % 
% end;
% 
% for i=1:ns;
% for j=1:nk;        
%     vcon2(i,j,:)  = vcon((i-1)*nk+j,:); 
%     p2(i,j,:)  = p((i-1)*nk+j,:);  
%     bp2(i,j,:) = bp((i-1)*nk+j,:);  
%     ibp2(i,j,:)= ibp((i-1)*nk+j,:);  
%     ct2(i,j,:) = ct((i-1)*nk+j,:);  
%     cn2(i,j,:) = cn((i-1)*nk+j,:);  
%     const2(i,j,:) = const((i-1)*nk+j,:);      
%     gn2(i,j,:) = gn((i-1)*nk+j,:);  
%     h2(i,j,:)  = h((i-1)*nk+j,:);  
%     m2(i,j,:)  = m((i-1)*nk+j,:);  
%     tau2(i,j,:)= tau((i-1)*nk+j,:);  
% 	rec2(i,j,:)= rec((i-1)*nk+j,:); 
%     def2(i,j,:)= def((i-1)*nk+j,:);  
% 	q2(i,j,:)  = q((i-1)*nk+j,:); 
% 
%     vaut2(i,j)  = vaut((i-1)*nk+j); 
%     paut2(i,j)  = paut((i-1)*nk+j);  
%     ctaut2(i,j) = ctaut((i-1)*nk+j);  
%     cnaut2(i,j) = cnaut((i-1)*nk+j);  
%     constaut2(i,j) = constaut((i-1)*nk+j);      
%     gnaut2(i,j) = gnaut((i-1)*nk+j);  
%     haut2(i,j)  = haut((i-1)*nk+j);  
%     maut2(i,j)  = maut((i-1)*nk+j);  
%     tauaut2(i,j)= tauaut((i-1)*nk+j);  
% end;
% end;
% 
% 
% w2    = alpha*p2.*h2.^(alpha-1);
% waut2 = alpha*paut2.*haut2.^(alpha-1);
% % 
% % sigma^2 = 0.029^2/(1-0.777^2)
% % mu = 0
% 

% 
% 
% if (nk==2);
% %default decision and prices
% def21(:,:) = def2(:,1,:);
% def22(:,:) = def2(:,2,:);
% 
% price21(:,:) = q2(:,1,:);
% price22(:,:) = q2(:,2,:);
% 
% figure;
% subplot(2,2,1);
% contourf(b(1:nb),a(1:ns),def21(:,:),20);
% xlabel('b');
% ylabel('a');
% title('default decision \kappa_{L}');
% subplot(2,2,2);
% contourf(b(1:nb),a(1:ns),def22(:,:),20);
% xlabel('b');
% ylabel('a');
% title('default decision \kappa_{H}');
% subplot(2,2,[3 4]);
% plot(b(1:nb),price21(lowa,:),'-r',b(1:nb),price21(meda,:),':r',b(1:nb),price21(hgha,:),'--r',...
%     b(1:nb),price22(lowa,:),'-b',b(1:nb),price22(meda,:),':b',b(1:nb),price22(hgha,:),'--b','LineWidth',1.5);
% xlabel('b');

% %%
% %value functions: repayment vs autarky
% vcon21(:,:) = vcon2(:,1,:);
% vcon22(:,:) = vcon2(:,2,:);
% vaut21(:,:) = vaut2(:,1);
% vaut22(:,:) = vaut2(:,2);
% 
% figure;
% subplot(2,1,1);
% plot(b(1:nb),vcon21(lowa,:),'-r',b(1:nb),vcon21(meda,:),':r',b(1:nb),...
%     vaut21(lowa,:)*ones(size(b)),'-b',b(1:nb),vaut21(meda)*ones(size(b)),':b','LineWidth',1.5);
% xlabel('b');title('value functions \kappa_{L}');
% axis([b(1) b(nb) -35 -15]);
% subplot(2,1,2);
% plot(b(1:nb),vcon22(lowa,:),'-r',b(1:nb),vcon22(meda,:),':r',b(1:nb),...
%     vaut22(lowa,:)*ones(size(b)),'-b',b(1:nb),vaut22(meda)*ones(size(b)),':b','LineWidth',1.5);
% xlabel('b');title('value functions \kappa_{H}');
% axis([b(1) b(nb) -35 -15]);
% end;
% %%
% %optimal policies and prices: repayment vs autarky
% NameArray = {'LineStyle','Color','LineWidth'};
% ValueArray = {'-',':','-',':';'red','red','blue','blue';1.5,1.5,1.5,1.5}';
% 
% NameArray0 = {'LineStyle','Color','LineWidth'};
% ValueArray0 = {'-','-';'red','blue';1.5,1.5}';
% 
% [ilowa_aut]=find(def2(lowa,1,:)==0);
% [ilowa_rep]=find(def2(lowa,1,:)==1);
% h2(lowa,1,ilowa_aut) = NaN;
% ha2(lowa,1,1:nb) = haut2(lowa,1);
% ha2(lowa,1,ilowa_rep) = NaN;
% 
% [ihgha_aut]=find(def2(hgha,1,:)==0);
% [ihgha_rep]=find(def2(hgha,1,:)==1);
% h2(hgha,1,ihgha_aut) = NaN;
% ha2(hgha,1,1:nb) = haut2(hgha,1);
% ha2(hgha,1,ihgha_rep) = NaN;
% 
% p2(lowa,1,ilowa_aut) = NaN;
% pa2(lowa,1,1:nb) = paut2(lowa,1);
% pa2(lowa,1,ilowa_rep) = NaN;
% 
% p2(hgha,1,ihgha_aut) = NaN;
% pa2(hgha,1,1:nb) = paut2(hgha,1);
% pa2(hgha,1,ihgha_rep) = NaN;
% 
% gn2(lowa,1,ilowa_aut) = NaN;
% gna2(lowa,1,1:nb) = gnaut2(lowa,1);
% gna2(lowa,1,ilowa_rep) = NaN;
% 
% gn2(hgha,1,ihgha_aut) = NaN;
% gna2(hgha,1,1:nb) = gnaut2(hgha,1);
% gna2(hgha,1,ihgha_rep) = NaN;
% 
% bp2(lowa,1,ilowa_aut) = NaN;
% bp2(hgha,1,ihgha_aut) = NaN;
% 
% cn2(lowa,1,ilowa_aut) = NaN;
% cna2(lowa,1,1:nb) = cnaut2(lowa,1);
% cna2(lowa,1,ilowa_rep) = NaN;
% 
% cn2(hgha,1,ihgha_aut) = NaN;
% cna2(hgha,1,1:nb) = cnaut2(hgha,1);
% cna2(hgha,1,ihgha_rep) = NaN;
% 
% ct2(lowa,1,ilowa_aut) = NaN;
% cta2(lowa,1,1:nb) = ctaut2(lowa,1);
% cta2(lowa,1,ilowa_rep) = NaN;
% 
% ct2(hgha,1,ihgha_aut) = NaN;
% cta2(hgha,1,1:nb) = ctaut2(hgha,1);
% cta2(hgha,1,ihgha_rep) = NaN;
% 
% m2(lowa,1,ilowa_aut) = NaN;
% ma2(lowa,1,1:nb) = maut2(lowa,1);
% ma2(lowa,1,ilowa_rep) = NaN;
% 
% m2(hgha,1,ihgha_aut) = NaN;
% ma2(hgha,1,1:nb) = maut2(hgha,1);
% ma2(hgha,1,ihgha_rep) = NaN;
% 
% w2(lowa,1,ilowa_aut) = NaN;
% wa2(lowa,1,1:nb) = waut2(lowa,1);
% wa2(lowa,1,ilowa_rep) = NaN;
% 
% w2(hgha,1,ihgha_aut) = NaN;
% wa2(hgha,1,1:nb) = waut2(hgha,1);
% wa2(hgha,1,ihgha_rep) = NaN;
% 
% tau2(lowa,1,ilowa_aut) = NaN;
% taua2(lowa,1,1:nb) = tauaut2(lowa,1);
% taua2(lowa,1,ilowa_rep) = NaN;
% 
% tau2(hgha,1,ihgha_aut) = NaN;
% taua2(hgha,1,1:nb) = tauaut2(hgha,1);
% taua2(hgha,1,ihgha_rep) = NaN;
% 
% figure;
% subplot(4,2,1);P = plot(b,squeeze(h2(lowa,1,:)),b,squeeze(ha2(lowa,1,:)),b,squeeze(h2(hgha,1,:)),b,squeeze(ha2(hgha,1,:)));
% title('h');xlim([b(1) b(end)]); ylim([0.7 1.01]);% -0.05 1.05]);
% set(P,NameArray,ValueArray);
% subplot(4,2,2);P = plot(b,squeeze(p2(lowa,1,:)),b,squeeze(pa2(lowa,1,:)),b,squeeze(p2(hgha,1,:)),b,squeeze(pa2(hgha,1,:)));
% title('p');xlim([b(1) b(end)]);%axis([b(1) b(end) 0.495 0.7]);
% set(P,NameArray,ValueArray);
% subplot(4,2,3);P = plot(b,squeeze(gn2(lowa,1,:)),b,squeeze(gna2(lowa,1,:)),b,squeeze(gn2(hgha,1,:)),b,squeeze(gna2(hgha,1,:)));
% title('g_{N}'); xlim([b(1) b(end)]);%axis([b(1) b(end) -0.05 1.05]);
% set(P,NameArray,ValueArray);
% subplot(4,2,4);P = plot(b,squeeze(bp2(lowa,1,:)),b,squeeze(bp2(hgha,1,:)));
% title('b''');xlim([b(1) b(end)]);%axis([b(1) b(end) -1.2 0]);
% set(P,NameArray0,ValueArray0);
% subplot(4,2,5);P = plot(b,squeeze(cn2(lowa,1,:)),b,squeeze(cna2(lowa,1,:)),b,squeeze(cn2(hgha,1,:)),b,squeeze(cna2(hgha,1,:)));
% xlim([b(1) b(end)]);%axis([b(1) b(end) 0.45 1.05]);
% set(P,NameArray,ValueArray);
% title('c_{N}');
% subplot(4,2,6);P = plot(b,squeeze(ct2(lowa,1,:)),b,squeeze(cta2(lowa,1,:)),b,squeeze(ct2(hgha,1,:)),b,squeeze(cta2(hgha,1,:)));
% title('c_{T}');xlim([b(1) b(end)]);%axis([b(1) b(end) 0.1 0.3]);
% set(P,NameArray,ValueArray);
% subplot(4,2,7);P = plot(b,squeeze(w2(lowa,1,:)),b,squeeze(wa2(lowa,1,:)),b,squeeze(w2(hgha,1,:)),b,squeeze(wa2(hgha,1,:)));
% title('w');xlim([b(1) b(end)]);%axis([b(1) b(end) -0. 0.3]);
% set(P,NameArray,ValueArray);
% subplot(4,2,8);P = plot(b,squeeze(tau2(lowa,1,:)),b,squeeze(taua2(lowa,1,:)),b,squeeze(tau2(hgha,1,:)),b,squeeze(taua2(hgha,1,:)));
% title('T'); xlim([b(1) b(end)]);%axis([b(1) b(end) -0.10 0.2]);
% set(P,NameArray,ValueArray);
% 
% %%
% % plot time series
% seriesp   = load('seriesp.txt');
% seriesr   = load('seriesr.txt');
% seriesh   = load('seriesh.txt');
% seriesb   = load('seriesb.txt');
% seriesct  = load('seriesct.txt');
% seriesm   = load('seriesm.txt');
% seriesgn  = load('seriesgn.txt');
% seriestau = load('seriestau.txt');    
% seriesa   = load('seriesa.txt');    
% seriesbor = load('seriesbor.txt');    
% 
% n = length(seriesp);
% 
% seriescn= seriesh.^alpha-seriesgn;
% seriesy = seriesa.*seriesm.^gamma+seriesp.*seriesh.^alpha; % -pm*seriesm;
% seriesgny = seriesgn./seriesy;
% 
% time =1:n;
% 
% ValueArray1 = {'-';'blue';1.5}';
% ValueArray2 = {'-','--';'red','blue';1.5,1.5}';
% ValueArray3 = {'-','--','-';'red','blue','black';1.5,1.5,1.5}';
% 
% % No shaded areas for default episodes
% % figure;
% % subplot(3,2,1); P=plot(time,seriesa');title('a');
% % set(P,NameArray,ValueArray1);
% % subplot(3,2,2);P=plot(time,seriesp);title('p');
% % set(P,NameArray,ValueArray1);
% % subplot(3,2,3);P=plot(time,seriesct);
% % %AX=legend('c_{T}','Location','NorthWest','Orientation','Horizontal');
% % title('c_T');
% % set(P,NameArray,ValueArray1);
% % %LEG = findobj(AX,'type','text');
% % %set(LEG,'FontSize',7);
% % subplot(3,2,4); PPP=plot(time,seriescn,time,seriesh,time,seriesgn);
% % axis([1 n -0.05 1.05]);
% % AX=legend('h','c_{N}','g_{N}','Location','NorthWest','Orientation','Horizontal');title('h,c_{N} & g_{N}');
% % set(PPP,NameArray,ValueArray3);
% % LEG = findobj(AX,'type','text');
% % set(LEG,'FontSize',7);
% % subplot(3,2,5); P=plot(time,seriesr);title('spreads');
% % set(P,NameArray,ValueArray1);axis([1 n -0.01 0.1]);
% % subplot(3,2,6); P=plot(time,-seriesb);title('b');
% % set(P,NameArray,ValueArray1);
% 
% 
% seriesbor = 1-seriesbor;
% seriesbor(199)=1;
% seriesbor(200)=0;
% seriesh(200) = 0;
% ntime = 200;
% 
% figure;
% subplot(3,2,1); 
% shadedTimeSeries(time,seriesa',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('a');
% AX1 = gca;box on;
% set(AX1,'YTick',0.85:0.1:1.15);
% axis(AX1, [time(1) time(ntime) 0.85 1.15]);
% subplot(3,2,2);
% shadedTimeSeries(time,seriesp',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('p');
% AX1 = gca;box on;
% set(AX1,'YTick',0:0.2:1);
% axis(AX1, [time(1) time(ntime) 0.95*min(seriesp) 1.05*max(seriesp)]);
% subplot(3,2,3);
% shadedTimeSeries(time,seriesct',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('c_T');
% AX1 = gca;box on;
% set(AX1,'YTick',0.05:0.05:1.15);
% axis(AX1, [time(1) time(ntime) 0.9*min(seriesct)  1.1*max(seriesct)]); 
% subplot(3,2,4); 
% shadedTimeSeries(time,seriesh',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('h,c_N & g_N');
% AX1 = gca;box on;
% axis(AX1, [time(1) time(ntime) 0 1.05]);
% hold on;
% H0 = line(time,seriesh','Color','b','LineWidth',2);
% H1 = line(time,seriescn','Color','r','LineWidth',2,'LineStyle','--');
% H2 = line(time,seriesgn','Color','g','LineWidth',2,'LineStyle',':');
% h = [H0, H1,H2];
% legend(h,'h','c_N','g_N','Location','SouthEast','Orientation','Horizontal');
% AX1 = gca;box on;
% axis(AX1, [time(1) time(ntime) 0 1.05]);
% subplot(3,2,5); 
% shadedTimeSeries(time,seriesr',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('spreads');
% AX1 = gca;box on;
% set(AX1,'YTick',0.0:0.025:0.1);
% axis(AX1, [time(1) time(ntime) 0 0.1]);
% subplot(3,2,6); P=plot(time,-seriesb);title('b');
% shadedTimeSeries(time,-seriesb',seriesbor',[],{[]}, [0.9 0.9 0.9], 3);title('b');
% AX1 = gca;box on;
% set(AX1,'YTick',0.0:0.02:ceil(max(-1.1*seriesb/0.001))*0.001);
% axis(AX1, [time(1) time(ntime) 0 ceil(max(-1.1*seriesb/0.001))*0.001]);